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We derive the Edgeworth streaming model (ESM) for the redshift space correlation function starting from 
an arbitrary distribution function for biased tracers of dark matter by considering its two-point statistics and 
show that it reduces to the Gaussian streaming model (GSM) when neglecting non-Gaussianities. We test 
the accuracy of the GSM and ESM independent of perturbation theory using the Horizon Run 2 A^-body halo 
catalog. While the monopole of the redshift space halo correlation function is well described by the GSM, 
higher multipoles improve upon including the leading order non-Gaussian correction in the ESM: the GSM 
quadrupole breaks down on scales below 30 Mpc/h whereas the ESM stays accurate to 2% within statistical 
errors down to 10 Mpc/h. To predict the scale dependent functions entering the streaming model we employ 
Convolution Lagrangian perturbation theory (CLPT) based on the dust model and local Lagrangian bias. Since 
dark matter halos carry an intrinsic length scale given by their Lagrangian radius, we extend CLPT to the coarse¬ 
grained dust model and consider two different smoothing approaches operating in Eulerian and Lagrangian 
space, respectively. The coarse-graining in Eulerian space features modified fluid dynamics different from dust 
while the coarse-graining in Lagrangian space is performed in the initial conditions with subsequent single 
streaming dust dynamics, implemented by smoothing the initial power spectrum in the spirit of the truncated 
Zel’dovich approximation. Finally, we compare the predictions of the different coarse-grained models for the 
streaming model ingredients to Al-body measurements and comment on the proper choice of both the tracer 
distribution function and the smoothing scale. Since the perturbative methods we considered are not yet accurate 
enough on small scales, the GSM is sufficient when applied to perturbation theory. 

PACS numbers: 


I. INTRODUCTION 

Redshift space distortions observed in galaxy surveys pro¬ 
vide a unique insight into the build-up of cosmological struc¬ 
ture by gravitational clustering of dark matter and its trac¬ 
ers such as halos and galaxies. Indeed, the redshift space 
two point correlation function carries valuable information on 
both, the real-space clustering and the peculiar velocity field 
since the observed redshift depends not only on distance but 
also on deviations from the overall Hubble flow. Peculiar ve¬ 
locities are generated by and hence correlated with the clus¬ 
tering of matter. 

There are two main effects in redshift space, a term intro¬ 
duced in [1], that affect the correlation function on large and 
small scales, respectively. On large scales the peculiar ve¬ 
locity associated with the coherent infall into overdense re¬ 
gions squashes structures and enhances the correlation func¬ 
tion along the line of sight which is captured by linear the¬ 
ory and known as the Kaiser effect [2]. On small scales, the 
elongation of nonlinear structures along the line of sight, the 
so-called ‘Fingers of God’ effect coined in [3] and first de¬ 
scribed in [4], leads to a suppression of the correlation func¬ 
tion. Based on this observation one of the first streaming 
models was developed in [5] by assuming an exponential rela¬ 
tive or pairwise velocity distribution with a scale-independent 
dispersion. Dispersion models [6, 7] aimed to phenomeno- 
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logically combine effects of linear clustering and small-scale 
velocity dispersion which act both multiplicative onto the 
redshift space power spectrum when their correlation is ne¬ 
glected. However, they have been shown to be unable to prop¬ 
erly account for redshift space distortions over a vast range of 
scales by means of 77-body simulations [8, 9]. 

To reunite the two disparate results for large and small 
scales, given by the linear theory [2] and the streaming model 
for nonlinear scales [5], the so-called Gaussian streaming 
model (GSM) was introduced in [10]. To obtain the GSM, 
the matter correlation function in redshift space was derived 
by considering the joint probability distribution of density and 
velocity. Assuming that the density is a Gaussian random held 
and the velocity is related to density as in linear perturba¬ 
tion theory one obtains a simple expression for the redshift 
space correlation function. It is given by a convolution of the 
real space correlation function and an approximately Gaussian 
pairwise velocity distribution whose mean and variance are 
given by the scale-dependent mean and variance of the pair¬ 
wise velocity. The GSM obtained via this approach can be un¬ 
derstood as generalization of the streaming model originally 
introduced in [5] to a scale-dependent rather than constant ve¬ 
locity dispersion which correctly reproduces the linear theory 
result [2]. The GSM, derived for the special case of Gaus¬ 
sian fluctuations in [10], has been generalized to fully non- 
Gaussian fields in [11]. Furthermore, therein a connection be¬ 
tween the redshift-space clustering and the pairwise velocity 
moments has been established. 

Furthermore, it has been shown recently in [12] that the 
assumption that the pairwise velocity distribution is locally 
Gaussian, with its mean and variance themselves Gaussian 
distributed allows to accurately recover the non-Gaussian 
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pairwise velocity distribution measured in simulations. This 
approach is different from the one presented here, where we 
assume that the mean and variance are not random variables 
but functions of separation that are either determined from 
data or inferred from theory. 

We start from a phase space distribution function for dark 
matter or its tracers, similar as done in [13, 14]. Indeed, our 
formulation relates the distribution function approach studied 
in Fourier space in [13, 14] and the Gaussian streaming model 
for redshift space distortions operating in configuration space. 
We decide to work in conhguration rather than Fourier space. 
A practical reason is that our formulation of perturbation the¬ 
ory will naturally produce expressions in real space. An¬ 
other argument is the fact that small spatial scales in the cor¬ 
relation function can be strongly affected by late-time bary- 
onic physics, while large scales, most importantly the batyon 
acoustic oscillations (BAG) peak, are not affected, see [15]. 
Therefore although late-time batyonic physics is confined to 
small r, it appears spread out in Fourier space. A generalized 
dispersion model, taking nonlinear couplings between density 
and velocity fields into account, has been proposed in [16] to 
provide consistent predictions for power spectra and correla¬ 
tion functions at the same time. 

To predict halo correlation functions in redshift space the 
GSM has been combined with perturbation theory to extract 
the streaming model ingredients, namely the real space cor¬ 
relation and the mean pairwise velocity and its dispersion, 
in [17-19]. A test of different analytic and phenomenologi¬ 
cal streaming models combined with perturbation theory, per¬ 
formed in [20], showed that they reasonably ht the simulations 
on intermediate scales 40Mpc//! < s < SOMpc/h while all 
models fail at small scales with Lagrangian schemes having 
the best performance around the scale of BAG. 

It is well known that no perturbative framework is able to 
accurately describe the fully nonlinear regime of structure for¬ 
mation. Fortunately, dark matter halos and their progenitors, 
which we denote by proto-halos, can themselves be treated 
as large cold dark matter (CDM) particles and therefore de¬ 
scribed by a pressureless dust fluid. The motion of these 
proto-halos is mostly determined by the large scale gravita¬ 
tional held and therefore much better describable with pertur¬ 
bation theory. The pressureless CDM Huid is described by a 
coupled system of differential equations consisting of conti¬ 
nuity, Euler and Poisson equations. These equations can be 
solved perturbatively - either in the Eulerian frame (SPT)[21] 
where everything is expanded in terms of density and velocity 
or in the Lagrangian frame (LPT) [22] where Huid-trajectories 
or displacement helds are considered. It is clear that the huid 
description should be applied only on scales larger than the 
particle size, in case of proto-halo “particles” this is the La¬ 
grangian size of the halos. Therefore it is natural to implement 
the Lagrangian halo size as a physically meaningful coarse- 
graining scale into the huid description for (proto-)halos [23]. 
This approach is to be seen in contrast to the so-called effec¬ 
tive held theory of LSS [24] for dark matter where the de¬ 
pendence of dark matter properties on the smoothing scale is 
unphysical and removed through renormalization. In order to 
model the trajectories of proto-halos we study in this paper a 


coarse-grained dust model in terms of the displacement held 
within Lagrangian perturbation theory (LPT). 

A big advantage of Lagrangian schemes [22] is the clearer 
physical picture they offer for the study of halo correlation 
functions, which are a key ingredient of the halo model [25] 
that is widely used in the analysis of galaxy, cluster and lens- 
ing surveys. In order to understand halo correlations one 
needs to understand the bias between the halo held and the 
underlying dark matter held. But halo bias is best under¬ 
stood using the spherical collapse model and excursion set 
theory [26, 27], both of which operate in the initial conditions 
and therefore in Lagrangian space, where they locally identify 
proto-halos within the initial density held and assign mass and 
collapse time to them. Therefore once the clustered or biased 
held of proto-halos is known it can be propagated to Eule¬ 
rian space using a Lagrangian method. Another advantage of 
Lagrangian methods concerns the convergence properties and 
the accuracy of the correlation function on the scales of in¬ 
terest, like the BAG scale or the mildly nonlinear scales. It 
is known that LPT performs much better on those scales, see 
the hrst Ligure of [28]; a higher precision is achieved with a 
smaller order in perturbation theory. The better convergence 
properties of the LPT displacement held compared to stan¬ 
dard perturbation theory (SPT) in Eulerian space are mainly 
due to fact that the relation between the density contrast and 
the displacement held is nonlinear and can be handled non- 
perturbatively. 

In hrst order LPT it is possible to analytically compute 
the density correlation function from the hrst order dis¬ 
placement held in a nonperturbative fashion which is called 
Zel’dovich approximation (ZA), see [29]. In the ZA par¬ 
ticles are displaced along straight trajectories, parametrized 
by the linear growth function, in a direction determined by 
their initial velocity. Despite its simplicity, the ZA is ca¬ 
pable of accurately describing gravitational dynamics over a 
surprisingly wide range of scales [28, 30]. In [30] the so- 
called truncated Zel’dovich approximation (TZA) was pro¬ 
posed as phenomenological method to improve the agreement 
between ZePdovich and proper A-body simulations by arti- 
hcially smoothing the initial power spectrum at the nonlin¬ 
ear scale of the hnal time of the simulation. The effect of 
the smoothing is to decrease the velocity in high density re¬ 
gions thereby reducing the amount of shell-crossing events 
and subsequent erasure of overdensities. Therefore, counter¬ 
intuitively, smoothing the initial power spectrum, which re¬ 
duces the initial power on small scales, actually can increase 
the hnal power on those scales. Locusing on statistical prop¬ 
erties of the nonlinearly evolved density held like the power 
spectrum, the TZA amounts to smoothing the linear initial 
power spectrum without affecting the dynamics itself. A de¬ 
tailed study and comparison between different hlters in [31] 
revealed that a Gaussian hltering scheme leads to best agree¬ 
ment with A-body data and considerable improvement over 
sharp k-truncation as originally suggested in [30] and top-hat 
hltering as studied in [32]. 

It is known that the Post-Zel’dovich approximation (PZA), 
where the displacement helds are calculated from second or¬ 
der LPT, improves over the ZA. Accordingly, the truncated 
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Post-Zel’dovich approximation (TPZA) with a smoothed ini¬ 
tial power spectrum performs even better than TZA, compare 
[33, 34]. We apply the framework of Convolution Lagrangian 
perturbation theory (CLPT) developed in [35] which recovers 
the ZA at lowest order while providing an approximation to 
PZA at higher order. CLPT can be understood as a partial 
resummation of the formalism presented in [36] providing a 
nonperturbative resummation of LPT that incorporates non¬ 
linear halo bias. We will compare two different smoothing ap¬ 
proaches within CLPT, namely a coarse-graining in Eulerian 
space (cgCLPT) with a coarse-graining in Lagrangian space 
implemented by smoothing the initial power spectrum in the 
spirit of the truncated Zel’dovich approximation (TCLPT). 
Those two procedures are distinct since a coarse-graining in 
Eulerian space also modihes the underlying dynamics becom¬ 
ing manifest beyond linear order in Lagrangian space, see 
[37], while our coarse-graining in Lagrangian space only af¬ 
fects the initial conditions. 

Structure This paper is organized as follows: In Sec. II we 
derive the Edgeworth streaming model (ESM) for the redshift 
space correlation function starting from an arbitrary distribu¬ 
tion function for biased tracers of dark matter by considering 
its two-point statistics and show that it reduces to the Gaussian 
streaming model (GSM) when neglecting non-Gaussianities 
in the pairwise velocity distribution. We then demonstrate the 
accuracy of the GSM and ESM on the basis of A-body simu¬ 
lations employing the Horizon Run 2 halo catalog. In Sec. Ill 
we built up on existing work and describe how the ingredi¬ 
ents of the streaming models can be inferred from the dust 
model and propose two different coarse-grained generaliza¬ 
tions of the fluid description. In Sec. IV we compute the real- 
space halo correlation function and the halo velocity statistics 
for the dust model employing Convolution Lagrangian Per¬ 
turbation Theory (CLPT) with two different coarse-graining 
schemes, an Eulerian (cgCLPT) and a Lagrangian (TCLPT) 
one. We conclude and describe possible further interesting 
lines of study in Sec.V. A list of abbreviations commonly 
used within this work can be found in App. A. 


II. EDGEWORTH STREAMING MODEL 

In order to infer predictions for the halo correlation func¬ 
tion in redshift space we use the Gaussian streaming model, 
originally derived in [10] and studied in [17] for the dust 
model. Starting from an arbitrary distribution of proto-halos 
we present a self-contained derivation of the Gaussian stream¬ 
ing model (GSM) from general assumptions which allows to 
include non-Gaussian corrections leading to the Edgeworth 
streaming model (ESM). We test the accuracy of the GSM and 
ESM using A-body simulation data from the Horizon Run 2 
(HR2) [38, 39] independent of perturbation theory. We then 
describe in Sec. Ill how the ingredients of the streaming mod¬ 
els can be inferred from the dust model and it’s coarse-grained 
generalization and present the CLPT computation and results 
in Sec. IV. 


A. Derivation of the ESM 

Let the phase space distribution function of dark mat¬ 
ter tracers X (like galaxies, clusters or halos) be given by 
fx(r, u, t). In this section we do not make any assumptions 
about its dynamics or statistical properties apart from that it is 
spatially statistically homogeneous 

{fx(ruUut)fx{r 2 ,U 2 ,t)) = {fxfx){r = r 2 - ri, iti, M 2 , f) ■ 

( 1 ) 

In addition we assume that the tracer density field 

1 -H 6x{r, 0^ J fx(r, u, f) (2) 

and higher moments Jd^u fx(r,u,t)ui^...Ui^^ are statistically 
homogeneous and isotropic. 

The observed position of a tracer Sobs - its angle on the sky 
ffobs and its observed redshift Zobs - corresponds to a point on 
the observer’s past light cone. As a first step towards calcu¬ 
lating tracer correlations on the past light cone, we will make 
two common simplifying assumptions. Eirst, since we are in¬ 
terested in equal-time correlation functions, we will approxi¬ 
mate the light cone in the neighbourhood of t by the f=const 
slice. Secondly, we use the distant observer approximation, 
where the line of sight is assumed to be a fixed direction z 
which is without loss of generality chosen as the direction of 
the z-axis, to relate the observed redshift-space position s of a 
dark matter tracer to its real-space position r . Those approxi¬ 
mations are despite their simplicity sufficient even for modern 
wide-area surveys within the level of current error bars, see 
e.g. Eig. 10 in [40]. Eor a general definition of redshift space 
and a discussion of wide-angle effects in linear perturbation 
theory we refer to [41]. In the distant observer approximation 
the observed comoving distance in redshift space s is affected 
by the peculiar velocity v ■ z - of the tracer along the line 
of sight via 

s - r+ ■ z) z , (3a) 

where 'Ll = aH = d and u - av. The observed position of the 
tracer perpendicular to the line of sight remains unaffected 
if we neglect gravitational lensing. In contrast, its coordinate 
i|l parallel to the line of sight z depends on the peculiar veloc¬ 
ity 

Sj_ - rj_ , i|| = s ■ i = ry H-. (3b) 

Since objects cannot disappear going from real space to red¬ 
shift space (assuming that all objects remain observable) we 
have the following relation between the densities in real and 
redshift space 

(1 + dx(«, 0) = (1 + ^x(r', 0) dV . (4) 

Although the correction to the real space position in redshift 
space is very small « ry, the clustering is affected con¬ 

siderably since the change of volume measure between real 
and redshift space, given by the Jacobian between d^s and 
d^r, involves the gradient of in linear perturbation theory 
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[2]. In the distant observer approximation, the tracer density 
fluctuation in redshift space (4) can be equivalently written as 


1 + 6 x(s 




d M fx(r, u, s - r - 


a^H 


( 5 ) 


which holds even for the case where the tracer velocity v is not 
a single valued function of r but instead has multiple streams 
or a continuous distribution, see also [42]. Later, in Sec. Ill, 
we will consider the special case of single streaming tracers 
described by the dust model for which this relation simplifies 
to (26a). 

We are interested in the redshift space two-point correlation 
function 


1 + 0 = (d + ^x(^i))(l + Sxisi))) , (6) 


where s — S 2 - si. By inserting (5) in (6) and re-expressing 
the delta functions in Fourier space and integrating over R - 
T] -H r 2 and one momentum variable the correlation function 
can be brought into the following form 


l+^x(s 




(2;r)3 ‘ 


Jk( 


J - (k ■ z) z,tj , 


Z(r, J,t) = j' d^Mi J" d' 


M 2 {fxfx}{r,Ui,U 2 ,t) 


X exp 


XU 2 - Ml) ■ J 


a^H 


(7a) 


(7b) 


where Z is the pairwise generating function. Next we Taylor 
expand W{J) := InZ around J = 0 


onto the line of sight k„ = ■ ■ -Zi,, (13). Expanding W 

in (9a) up to second order in J implies that all redshift space 
distortion induced clustering is encoded in the scale dependent 
mean and variance given by the pairwise velocity v \2 and its 
dispersion cTj^. As we will shortly see, this corresponds to the 
Gaussian streaming model (GSM). Since the GSM is known 
to be a good approximation, we will perform an expansion 
around this Gaussian 


ln(l + ^xir, t)) + ivx 2 ■ ^ - 2 ^ 


exp 

~ (1 +^x(Lf))exp 


iV\2 ■ J - - J CTjjJ 


1 - 
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This approach is similar to the idea behind Convolution La- 
grangian perturbation theory (CLPT), see [35]. To obtain the 
Gaussian streaming model it is crucial to expand in cumu- 
lants and keep the pairwise velocity mean and dispersion in 
the exponent, corresponding to specific resummation of mo¬ 
ments. Within the distribution function approach to redshift 
space distortions developed in [13, 14, 42] a moment expan¬ 
sion without such an resummation was performed such that 
the connection to the Gaussian streaming model is not mani¬ 
fest and has not been discussed. 

Later when testing the accuracy of this model, we will re¬ 
strict ourselves to the leading order non-Gaussian term. How¬ 
ever one can systematically expand the exponential of the non- 
Gaussian contributions to Z(J) in an Edgeworth expansion 
[43, 44] around a Gaussian pairwise velocity probability dis¬ 
tribution. The Edgeworth series E„ is an asymptotic expansion 
to approximate a probability distribution using its cumulants 
K„. With the Gaussian distribution as reference function it can 
be written as, see Eq. (43) in [45], 


OO J 

W(J) = y , Kn 

n\ 


d"W 


n=0 


{dijy 


( 8 ) 


y=o 


Keeping only the terms up to third order n = 3 we obtain 


En{x) = 


xI27TK2 


exp 


{X-Kif 


‘-ZZ 


2 k 2 

■■■As-r+A 


(11a) 


Rs+lr 


s=l r=\ 


X - K\ 


VT( J) ^ ln(l -H ^x{r, t)) + ivn ■ J ■ 
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with the cumulants k„ as expansion coefficients 


1 




(9a) where A„ are the normalized and rescaled cumulants 

Kfj 


(9b) 


Afj — 


n{n - l)/ty 

the Bell polynomials 


(11b) 


1 H- ^x(r, t) exp A-o = Z |y=o . 

dz I 

(5/J) l/=0 


v\ 2 (r,t) /«i = 


(r,f) K 2 = 


(l+^x(Lf)) 

d^Z I 
{idJ)- i/=0 


dZ dZ 
{idJ) {idJ) 1 7=0 


(1+^x(l0) {l+^x{r,t)A 


= d-j2(r, f) - vn{r, f)vx2(r, t ), 


A;f(r,f):=xf = 


12 


(0-2 


(10a) 

(10b) 


(10c) 

(lOd) 


whereAl'lB*' +A’^‘BX Since we have to evalu¬ 

ate all expressions atJ- (k-z) z we project the cumulants k„ 


Til, 1 ( 43 ) - /I 3 , By, 1 ( 43 ,/I 4 ) - 44 , £ 2 , 2 ( 43 ) - 43 , ( 11 c) 

and Hn the probabilists’ Hermite polynomials 

EIAx) = x^ -3x , HAx) - x^ - 6x2 ^ ^ 

H(,(x) = X® - 15x'^ H- 45x2 - 15 . 

In the following we perform the Edgeworth expansion up to 
n = I explicitly taking into account the first non-Gaussian 
correction given by the pairwise velocity skewness = A 12 . 
Kurtosis K 4 would arise in the next order but won’t be consid¬ 
ered in this paper. 
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We can now plug the Edgeworth expansion (11a) of 
Z = exp W according to Eq. (9) into the correlation function 
Eq. (7a). In the course of the calculation we will use cylindri¬ 
cal coordinates 


s — sx[cos(0)& + sin(0)y] + s\\z , 


since ^x{s, t) does not depend on the angle <p. Performing five 
of the six integrals in Eq. (7a) we obtain the Gaussian stream¬ 
ing model (GSM) Eq. (12a) at second order in the cumulant 
expansion and the leading and up to order n corrections of the 
Edgeworth streaming model (ESM) (12b), (12c) at third and 
n-th order, respectively 


1 

X exp 


f 


drii 


-oo V^cri 2 (r, ry, f) 

(S|l - ry - Vi2(r,t)rii/rf 




Ai2 
1 + — 




/ ^srv \ 
1^12 ' 


n=oo s 

zz 

.9=2 r=l 




H, 


s+ 2 r 


^srv 
O' 12 


(12a) 

(12b) 

(12c) 


In more detail, the integral in Eq. (7a) introduces the pair¬ 
wise probability distribution multiplied by (1 +^x(r, t)), while 
the trivial k;t,ky integrals enforce The integral 

ensures while the 0-integral gives a factor of Itt. 

We defined 


vnir, t)r\\/r t!i2(r, t) ■ z ^ ki , 

A.,™ sy - ry - t;i 2 (r,f)r||/r, 
cri 2 (r, ry, f) i^crj 2 (r, t)z = ki 

= <5-i2(r, ry, f) - yi 2 (r, tf{r\\lrf 
= cr2(r, t){r\\lrf + o-^(r, t) [l - (ry/r)^] , 
Ai2 = K2 

= (A||(ry/r)^ + A^ [l - (ry/r)^]) ry/r. 


(13a) 


(13b) 

(13c) 

(13d) 


In a previous study of the GSM [17], the following formula, 
inspired by the exact result from [10] for the case where both 
density and velocity fields are Gaussian and related to one an¬ 
other as in linear theory, was suggested to calculate Gaussian 
streaming redshift space distortions 


1 -Hfz(sii,s±,f) = r 

v/—cx 


X exp 


-(l+a(r,f)) 

V27ro-i2(r, ry, f) 

(■^11 ~ Ol ~ Vi2ir,t)r\\lrf 


where our sy and ry corresponds to ry and y used in [17, 18], 
respectively. Note that, (14) corresponds to (12a) when the 
variance, given by the second pairwise velocity moment crj 2 , 
is replaced by the pairwise velocity dispersion The two 
quantities are related via Eq. (13c) such that cr^ = 


and cr^ = d\. By expanding Z ^ exp(Wo + Wi + W 2 ) one 
obtains the GSM (12a) with the second cumulant cr ^2 vari¬ 
ance whereas when expanding Z ^ Zo(l -H Zi/Zo + Z 2 /Z 0 ) - 
Zo exp(Zi /Zo + Z 2 /Z 0 ) one obtains the GSM (14) with the sec¬ 
ond moment (Xjj as variance. When linearized, both expres¬ 
sions (12a) and (14) agree, because iij 2 second order, and 
correctly reproduce the Kaiser formula as shown in [10, 17]. 



FIG. 1: Comparison between the Horizon Run 2 measurement {data 
points) of the multipoles of the redshift space correlation function de¬ 
fined in (15) and the GSM (12a) using {thick dashed) or {thin 
dashed) related via (13c). upper panel The quadrupole ^2 shifted by 
10(i - l)(Mpc//!)^ where i labels the mass bin. middle panel The 
hexadecapole ^4 shifted hy 100(i - l)(Mpc//!)^. lower panel The 
hexacontatetrapole (64-pole) ^6 shifted hy 20(i - l)(Mpc//!)^. 

It is natural to follow an expansion in W and to keep only 
the Gaussian part in the exponential in case the pairwise ve¬ 
locity distribution is close to a Gaussian. On the other hand 
the moment expansion of [42] is natural from a perturbation 
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theory perspective, in which only moment spectra are kept that 
are nonzero up to certain order in perturbation theory. 

That the pairwise distribution function is indeed approxi¬ 
mately Gaussian with a variance given by (T 12 rather than 0-12 
becomes clear in Fig. 1, where we compare the GSM with 
the second cumulant 0-^2 (12a) to the GSM with the second 
moment cfj 2 (14) as the variance of the Gaussian. The exact 
definition of the redshift space multipoles depicted in Fig. 1 
and the reason for their normalization will be given in the next 
subsection. As we can clearly see the use of the second cumu¬ 
lant (Tj 2 significantly improves the agreement for the redshift 
space distribution function with the A-body simulation com¬ 
pared to the model where the second moment crj 2 is used. In 
[17] it has been phenomenologically accounted for that differ¬ 
ence by subtracting the square of the mean infall u ^2 from 
to get the dispersion about the mean. We leave it for future 
work to directly compare the ESM to the distribution function 
approach [14]. 


B. Accuracy of the GSM and ESM 

In the following, we assess the accuracy of the GSM (12a) 
and the leading order of the ESM (12b) by comparing the re¬ 
sults of the corresponding integrals ( 12 ) with the directly mea¬ 
sured redshift space halo correlation function f(s,ju, f). This is 
done by inserting the real space correlation ^(r), the pairwise 
velocity vnCr) and velocity dispersion crj 2 (r,yu) measured in 
an A-body simulation into Eq. (12a) and additionally measur¬ 
ing the skewness Ai 2 (r,/i) and plugging it into Eq. (12b). 

The Horizon Run 2 (HR2) A-body simulation [38, 39] has 
an enormous size of 7200Mpc//! and consists of 6000^ parti¬ 
cles of mass IgM = 11.097. Eor the mass units we use the 
notation IgM s log[o(M/!/Mo). We measured halo correla¬ 
tion functions and velocity statistics from large galaxy-sized 
haloes IgM = 13.0 to cluster-sized halos IgM = 15.2 at the 
redshift z = 0. In an accompanying work [46] we describe 
in detail how the correlation functions and Gaussian stream¬ 
ing ingredients haven been determined from the HR2 halo 
catalog. In order to evaluate and compare ^x('S||j Sr- 0 from 
Eq. (12) to simulations it is useful to expand fx('S||,'S±) into 
Legendre polynomials L„(/2) using -H and n - sy/s 

CO 

^x{s,li) = ^ L„(^)^x,n(s ), (15a) 

n =0 

^x(s,fi,t)Ln(jJ.)dii. (15b) 

vanishes for all odd n. In linear perturbation theory, 
the only non-zero multipoles are the monopole ^0, quadmpole 
^2 and hexadecapole ^4 and even in the nonlinear regime the 
magnitude of rapidly decreases with n. The linear results 
go back to [47] and are given in [17] for the case of Eulerian 


bias as 

^ois) = + ^b^f + ^ J dkk^PLik)jo{ks) 

^^(s) = - ^^b^f + ^ Jdk^P^miks) (16) 

^^(■s) = dkk^PL(k)j4(ks), 

where / is the linear growth rate and j„(x) are the spherical 
Bessel functions. We use their prefactors, given in terms of 
linear local Eulerian bias b^ - I H-7>i(lgMQpj) determined from 
the best fitting mass for the real space correlation function, 
see Tab. I, as a normalization when plotting multipoles. In 
Eig. 2 we compare the redshift space halo correlation function 
predicted from the GSM (12a) and ESM (12b), by measur¬ 
ing their ingredients from the HR2 data, to the direct mea¬ 
surements within HR2 for the redshift-space monopole ^ 0 , 
quadrupole ^2 and hexadecapole ^ 4 . We find that the ESM 
( 12 b) clearly improves the quadrupole ^2 and hexadecapole ^4 
on small scales compared to the GSM (12a). As evident from 
Eig. 3 the quadrupole predicted by ESM (12b) is accurate to 
2 % within statistical errors down to lOMpc//; in contrast to 
the GSM (12a) which breaks down below 30Mpc//!. A simi¬ 
lar trend can be observed for the hexadecapole ^4 whereas the 
monopole ^0 is less sensitive to non-Gaussian terms. Appar¬ 
ently smaller halos are more sensitive to non-Gaussian correc¬ 
tions which is in line with the expectation that smaller objects 
are more affected by nonlinear dynamics. 

We conclude that the GSM is a very accurate model for the 
multipoles n = 0 , 2 ,4 of the redshift space halo correlation 
function on scales larger than 30Mpc//! while the ESM stays 
accurate down to lOMpc/Zz. Our result is consistent with the 
previous finding that the GSM monopole is accurate on the 
percent level down to lOMpc//? and the quadrupole down to 
30Mpc//!, compare Eig. 6 in [17]. This shows that the expan¬ 
sion of Z around J - 0 was justified and that halos over a 
wide range of masses can indeed be reasonably described by 
the GSM/ESM (12a/12b). 

Having established the range of validity of the streaming 
models GSM and ESM, we can use them as a basis for the 
theoretical modeling of redshift space halo correlation func¬ 
tions being aware of their limitations. As a next step, accurate 
theoretical predictions for the streaming model ingredients, 
i>i 2 (r), crj 2 (r,ju) and Ai 2 (r,ju) are needed. In the follow¬ 
ing two sections we will combine the GSM/ESM with pertur¬ 
bation theory employing that halos can be treated as single¬ 
streaming objects when the fluid description is only applied 
on scales larger than their size, given by the Lagrangian ra¬ 
dius. More precisely, we will calculate the streaming model 
ingredients from Convolution Lagrangian Perturbation The¬ 
ory (CLPT) based on the dust model and extend it to include 
a coarse-graining scale chosen to be the Lagrangian radius. 

C. Pairwise generating and tracer cumulants 

By performing a cumulant expansion we can relate the term 
J ■ (U 2 - Ml) contained in the exponential of the generating 



7 







FIG. 2: The redshift-space multipoles ^2 and ^4 at z = 0 predicted by the GSM (12a) (thin dashed) and ESM (12b) (thick dashed) using 
the HR2 data compared to the direct HR2 measurement (data points) normalized with respect to their bias factors for the different mass bins. 
Similar to previous plots we added a mass bin dependent constant to all curves for better visibility, upper left panel The monopole ^o(‘S) shifted 
by 10(i - 1)( Mpc//;)-. upper middle and right panel The quadrupole ^ 2 (^) shifted by 20(i - 1)( Mpc/h)^. lower panel The hexadecapole ^ 4 ( 5 ') 
shifted by 200(i - l)(Mpc//i)^ 


function Z from (7b) to the cumulants of the tracer distri¬ 
bution function fx- Therefore we introduce the moment gen¬ 
erating functional 


G[J] 



exp 


\iJ ■ It] fx{r, u, t), 


(17) 


(ti) 

which allows to compute the cumulants Cy of the distribution 
function fx according to 


This can be used to re-express Z (r, J) from (7b) 


Z = 


fx,lfx,2 exp [ij ■ (U 2 - Ml)] 

fx,i exp [-iJ ■ Ml] 

- (G[J](r 2 )G[-J](ri)) (19) 

“ ;N j. j. 

V I Jii-JiN (ym ,/I'vA'ptA') ) 

Z_i „2NjjN V^X,2,ii-iN ' ^ 


p'in) 

X,ivi„ 


:= (-0‘ 


5"lnG[J] 

^Ji^ ... dJi^^ 


(18) 


where J = JI(a^H), fx,i( 2 ) = /x(r'i(2), Mi(2), 0 and C^xm 2 ) = 

Cf\rn2)). 


J=0 
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FIG. 3: The ratio of the GSM (12a) {lower panel) and ESM (12b) {upperpanel) prediction using the HR2 data and the direct HR2 measurement 
{data points) for the lowest three redshift-space multipoles, left panel The monopole fo{s) middle panel The quadrupole ^ 2 (''')- right panel The 
hexadecapole ^ 4 (i). 


III. DETERMINING STREAMING MODEL 
INGREDIENTS BASED ON THE DUST MODEL 

In this Section we describe how the scale-dependent func¬ 
tions entering the streaming model can be determined once a 
phase-space distribution fx of the tracers is specified. To draw 
conclusions based on theoretical modeling it is due to connect 
the (proto-)halo distribution fx to an underlying dark matter 
distribution / whose dynamics is known to be governed by 
the Vlasov-Poisson equation. Halos are biased tracers of dark 
matter, since according to spherical collapse and excursion set 
theory [26, 48, 49], the probability of forming a halo depends 
on the initial density field. Therefore, there are two steps for 
determining streaming model ingredients; 

1. Choose a model for the distribution function / of dark 
matter that reasonably approximates Vlasov dynamics. 

2. Specify a bias model in order to relate the (proto-)halo 
cumulants to the ones of dark matter 

In the following we employ the pressureless fluid model as 
standard model for cold dark matter and discuss different pos¬ 
sibilities to incorporate a coarse-graining in this fluid picture. 
For relating the halo to the dark matter density we use local 
Lagrangian bias with zero velocity bias and present two pos¬ 
sibilities to generalize this notion to higher cumulants. 


A. The single-stream case: dust model 

In the context of analytical modeling COM dynamics, usu¬ 
ally the dark matter distribution is assumed to be described by 
the pressureless fluid (dust) model 

/d(r, u,t) - {\ + 6{r, t))5Y){u - av{r, t)) (20) 

which encodes all properties in terms of a number density 
n{r) = 1 H- 6 {r) and a single-streaming and curl-free veloc¬ 
ity v{r) fulfilling the coupled continuity, Euler and Poisson 


equations [50]. The cumulants of the dust model are 

= In (1 -H d) , Cf'’ = avi , = 0 (21) 

which displays that the dust model is entirely described by 
density and velocity and all higher cumulants such as velocity 
dispersion vanish identically. Although the dust model is an 
exact solution of the Vlasov equation, its applicability is lim¬ 
ited to the single-stream regime. It does not allow to describe 
the nonlinear stage of structure formation during which higher 
cumulants are sourced by the occurence of shell-crossing after 
which multiple streams form. For (proto-)halos this limitation 
is not as severe since they approximately behave as single¬ 
streaming objects even though a large fraction of dark matter 
particles resides in halos where it is multi-streaming and not 
accessible by the dust model. Hence, the proto-halos can also 
be described in terms of a single-streaming dust fluid 

fxir, M, f) = (1 + Sx(r, t))dD (« - avx(r, t)) . (22) 

To connect the density of halos to the dark matter density, we 
assume local Lagrangian bias 

(1 + 6 x(r, t)) dV = Fmq), t] d^q . (23) 

This equality states that proto-halos identified in the linear ini¬ 
tial conditions, depending only on the smoothed initial linear 
density held 6 R{q), are conserved until they form a proper halo 
at time t. The proto-halo initial density held is assumed to be 
a local function F[6R{q), f] of the initial linear density held 
6 L{q) smoothed over some scale related to the Lagrangian 
size R of the proto-halo by applying a window function W 
in Fourier space 

/ j3 7 

-^W{kR)e‘^-^6L{k). (24) 

The choice of the appropriate smoothing scale R will be elabo¬ 
rated in more detail in a forthcoming paper [46]. Note that for 
the computations in both iPT and CLPT this smoothing scale 
R is effectively removed by setting the window function to 
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unity. In [36] this is justified by claiming that the large-scale 
clustering of biased objects should not depend on the artificial 
choice of R to define the background field and seconded by 
the assertion that this is demanded by consistency with the ap¬ 
proximation being valid only on scales larger than the smooth¬ 
ing radius R. In [18] it is furthermore argued that R naturally 
drops out in the final statistics of interest and is only necessary 
to keep intermediate quantities well-behaved. We will pre¬ 
serve the smoothing and see in Sec. IV, in particular Figs. 4 
and 5 how large the effect on the baryon acoustic peak is when 
a smoothing at the Lagrangian scale is performed compared to 
the case where the smoothing is dropped. 

The tracer velocity field Vx — a\['x displaces the proto¬ 
halos to their halo virialization sites r - q + ^x{q)- We 
assume zero velocity bias such that proto-halos move along 
dust fluid trajectories r = q + ^(q) with the dust velocity 

Vx - V - 


CLPT ; fxir, it, f) = (1 + dx(r, fjjdo (« - av{r, t)) 


(25a) 


Local Lagrangian bias (23) allows us to relate the densities in 
real and Lagrangian space in the following way 

1 H- 6x(r, 0 =J' F[dR(q), t]6D(r-q- ^(q, t)) . (25b) 

Due to the single-valuedness of the velocity we recover from 
inserting (25a) into (4) a simpler relation between densities in 
real and redshift space 


M) = |. 


1 -H (5x(s, f) = dv (1 -H 6x{r, 0) 5d U - r - 


v{r, t) 


aH 


(26a) 


We could combine both relations in a single expression by 
expressing the velocity in terms of the displacement v - a\P 


1 -H dx(s 




d q F[ 6 R(q),t]x 


(26b) 


X do I s - q - ^{q, t) - 


Note however that we will not use formulas (26) explicitly. 
Instead we will rely on the GSM (12) to go from real space to 
redshift space and Eq. (28) below, to go from Lagrangian to 
Eulerian space. 

Derivation We already showed in (12a) that one can ob¬ 
tain the GSM (14) from quite general assumptions, in particu¬ 
lar that no assumptions about tracer dynamics and bias are re¬ 
quired. Now, we will specialize Z from (19) to the dust ansatz 
(25) for the tracer phase-space distribution fx combined with 
local Lagrangian bias (23) as considered in [18]. 

Eirst, we use the dust model cumulants (21) applied to the 
proto-halo distribution from CLPT (25) and plug them into 
the general expression for Z(r, J) in terms of tracer cumulants 
(19) obtaining 

Z = / [1 H- 5x(r’i)] [1 + Sx(r 2 )] exp 


J ■ (v(r 2 ) - v(ri)) 


aH 


We then switch to Lagrangian space making use of local La¬ 
grangian bias (25b) and express the bias function 

^ (27) 

as well as the delta function 6v>(r - q - ^{q, t)) in Eourier 
space. Next we replace the single streaming velocity by the 
derivative of the displacement field v{r) - a^{q) and inte¬ 
grate over Q = qi + q 2 to obtain 


Z(r, J, f) 



ik-(q-r) rddi_d^ 

(27Tp J (iTTp 


X F(Ai)F(A2)(e‘^), 


(28a) 


with 


X — Ai 6 i + A 2 S 2 + k • ^ + J • 


A 

77 ’ 


(28b) 


where q ^ q 2 - qi, 6 ^ 2 ) = dR(qi( 2 )) and A = ^(q 2 ,t) - 
^{qi,t)- 

Previous studies Originally, in [35], Eq. (26b) was used to 
derive an expression for the two-point correlation function 


l+a(s,t) 


/•'’“I 

f 




{ 2 np 
d/lj d/l2 
( 2 np 


ik-{q-s) 


(29) 


■F{Ax)F{A2){e'^^-^-^’"'^^ , 


which was then evaluated within CLPT to obtain a Post- 
Zel’dovich approximation for biased tracers in redshift space. 
The formula (14) was suggested in [17] to calculate Gaussian 
streaming redshift space distortions, following the idea of [10] 
to reconcile the streaming model [5] for nonlinear scales with 
linear theory [2] by considering a scale-dependent variance. 
In [17], the pairwise velocity mean v \2 and second moment 
(Tj 2 entering the streaming model (14) were calculated from 
SPT with linear bias while the real space correlation ^(r) was 
inferred from LPT with local Lagrangian bias (25b). Later on, 
in [18], the real space correlation and velocity statistics were 
treated on the same footing and determined within CLPT [35] 
together with local Lagrangian bias. Note that (29) involves 
a three-dimensional q-integral which needs to be evaluated 
numerically within CLPT [35]. Studying the expression (28) 
for Z in CLPT, the streaming model ingredients can be calcu¬ 
lated according to (10) and involve at most two-dimensional 
numerical integrals [18]. This a practical reason to chose to 
perform an Edgeworth expansion of Z to obtain the Gaussian 
streaming model (12a) and its non-gaussian generalization - 
the Edgeworth streaming model - whose numerical evaluation 
is more efficient than the full CLPT expression (29). 


B. Beyond single-stream: coarse-graining the dnst model 

In the following we compare several distinct approaches 
of coarse-graining a dust fluid, namely a coarse-graining 
in Eulerian space (cgCLPT) and a coarse-graining in La¬ 
grangian space implemented by smoothing the initial power 
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spectrum in the spirit of the truncated Zeldovich approxima¬ 
tion (TCLPT). A key question is how to generalize the biasing 
scheme employed for CLPT based on the dust model to the 
coarse-grained case. So far we assumed local Lagrangian bias 
for the density (25b) and zero velocity bias. This might be 
generalized by (a) assuming zero velocity bias and that higher 
cumulants for the tracer vanish identically s 0 moti¬ 

vated by the fact that proto-halos can be described well by 
single-stream physics such that, in analogy to the CLPT case, 


fx(u, r, f) = (1 + 6x{r, t)) 6u(u - av(r, t)). (30a) 

or (b) assuming that tracers and dark matter are only biased 
with respect to density such that all higher tracer cumulants 
are identical to those of dark matter = ^(^'21) and 

. (30b) 

l+ 6 x(T,t) 

Note that in order to write the biasing in analogy to the 
dust case (28a) it is necessary that fxiu, r,t)l{l + 6x{r, t)) is 
independent of the bias function F which is achieved by both 
relations (30). Then, the redshift space correlation takes the 
form 


/Sr'"’"/ 


(27r)3 
X F{Ai)F(A2){e‘^}, 


d/lj d/l2 


(31a) 


with 


^ — Tidi -i- A 2 S 2 + ^ 


(31b) 


“ ;N-l J J 

N=1 


c 


m 




where C^xi{ 2 ) ~ (’'i(2)(9i(2)))- Hence, the ESM ingredi¬ 

ents are still computed according to Eqs. (10) with Z from 
(28) replaced by Z from (31). If we consider the GSM, ex¬ 
panding up to second order in J, we see that the first cumu- 
lant corresponds to the term A./H that is also present in 
the single streaming Gaussian streaming model (28b) and con¬ 
tributes both to the mean vn and variance cr ^2 of the Gaussian. 
In contrast, the second cumulant Cy is conceptually new and 
contributes only to the variance of the Gaussian, whereas all 
higher cumulants are irrelevant for the GSM but only 

contribute to the ESM. 


7. Coarse-graining in Eulerian space (cgCLPT) 

Coarse-graining the dust model on a length scale o-^ in Eu¬ 
lerian space and a velocity scale 01 , gives rise to the so-called 
coarse-grained dust model as described in detail in [37]. We 
shortly recap the main results that are of direct relevance here. 
The coarse-grained dust model is defined as a smoothing of 


the dust phase space distribution with a Gaussian filter of 
width cPc and 01 , in x and u space, respectively 


/d = 


/ 

f 


d^x d^M 
(27roicr„)3 

d^x 

(27r(7icr„)3 


exp 


exp 


(x - x)^ 

2 cry 

(x - x)^ 

2 cry 


(u - u)^ 

2cr„2 


fd(x, u) 

(32) 


(u - avixyf' 

2cr„2 


n{x). 


If .Ttyp and Mtyp are the (minimal) scales of interest we have 
to ensure that ctx JCtyp and ai, « Wtyp in order to be able to 
resolve these scales. The coarse-grained dust model features 
higher cumulants which are absent in the pressureless fluid 
case and given by 


C<‘’^ = lnh, 


nvj 

h 


avi , 


/^(2) _ 2 r 

C- ■ - oil ^ij ' 


ViVj nvi nuj 


cS = 


n 

•fcyc. perm. 


nUiVjVic nViVifiVk nvinvjfivk 
- - - -H2-T - 


(33a) 

(33b) 

(33c) 


see [37], where we also dehned 
d^x 


g(x) 


-f 




exp 


(x — 5)2 
2c7i2 


g(x). 


The coarse-grained velocity v is the mass-weighted dust ve¬ 
locity which is obtained by smoothing the momentum field 
nui and then dividing by the smoothed density field «. From 
a physical point of view v describes the center-of-mass ve¬ 
locity of the collection of particles inside a coarsening cell 
of diameter around x. Whereas the contribution from the 
velocity smoothing scale has no dynamical effect but only 
contributes to the velocity dispersion, the smoothing on a hxed 
Eulerian scale 01 affects the fluid dynamics [51,52]. The mod¬ 
ified fluid equations resulting from (32) and the dependence of 
the smoothing scale o} have been studied perturbatively in Eu¬ 
lerian and Lagrangian space [37]. The displacement field 
was defined as the integral lines of v =; and can be deter¬ 
mined perturbatively from the coarse-grained Eulerian quanti¬ 
ties S and V, as described in [37]. Our Eulerian-coarse-grained 
displacement ^ should not be confused with the direct coarse- 
graining of ^ in Lagrangian space considered in [24]. 

If we combine the biasing scheme (30a) for single¬ 
streaming tracers with the coarse-grained dust model (scg- 
CLPT), we effectively consider 


scgCLPT : fx(r, t6, f) = (l + Sx) (u - av(r)) 


(34a) 


which is analogous to (25) but the velocity and tracer density 
are now expressed in terms of the coarse-grained displacement 
field, V - and 

1 H- 6 xir) = J" d^q F[6 r=o-M)] Su{r-q- ^(q)) ■ (34b) 
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Then, the evaluation of the streaming ingredients (10) from 
(31) within cgCLPT is affected by the coarse-graining only via 
modified expressions for the kernels involving coarse-grained 
displacements We present the Lagrangian correlators for 
the coarse-grained dust model that are relevant for the CLPT 
evaluation in App. B. 

When instead choosing the biasing scheme (30b), we as¬ 
sume that halos and dark matter are only biased with respect 
to density 


cgCLPT : fx(r, u, t) = ^ f /^(r, u, t) 

1 -H 6 (r, t) 


(35) 


Then, in addition to having modified displacement kernels, 
the computation of the streaming ingredients is also affected 
by the occurrence of velocity dispersion encoded in higher 
cumulants (33). This affects the variance through 
as computed in App. C and the leading-order non-Gaussian 

-( 2 \ "CS) 

correction A^ through Cy and . 


2. Coarse-graining in Lagrangian space (TCLPT) 


Our coarse-graining in Lagrangian space is based on a fixed 
smoothing scale in the Lagrangian or initial condition space 
q corresponding to the scgCLPT model (34) 

TCLPT ; fxir, m, f) = (l + 5x,o-,) (m - aVo-^ir)) , (36a) 

but with the velocity and density given by the smoothed dis¬ 
placement field via Dq-, = 

1 + Sx,a-,{r) = d^q F[d/i=^^(q)] 5d (^ - 9 - ^o-,(g)) ■ (36b) 

If one considers coarse-grained dust dynamics in Lagrangian 
space then one is lead to [24]. For TCLPT, we instead per¬ 
form a smoothing in the initial conditions while keeping the 
pressureless fluid dynamics unchanged. This is implemented 
by smoothing the initial linear density field diiq), or equiv¬ 
alently the power spectrum Pi —» Pi when calculating the 
Lagrangian correlators for CLPT (25) that are given in [18] 


TCLPT; CLPT with Pi{k) exp Pi{k) 


Note that when evaluated in first order LPT in which case 
CLPT is identical to the Zel’dovich approximation), then also 
TCLPT (37) and cgCLPT (34/35) are identical to the TZA and 
differences arise when nonlinearities in the displacement field 
are taken into account. In the following section we evaluate 
the streaming model ingredients within CLPT up to second 
order in the power spectrum. 


IV. EVALUATION OF STREAMING MODEL 
INGREDIENTS WITHIN CONVOLUTION LAGRANGIAN 
PERTURBATION THEORY 


In this Section we evaluate the scale-dependent functions 
entering the streaming model within Convolution Lagrangian 


perturbation theory (CLPT), introduced in [35] based on the 
pressureless fluid model (25) and its coarse-grained versions 
scgCLPT (34), cgCLPT (35) and TCLPT (37) defined in the 
last section. First, we present the calculation of the real space 
correlation function and the pairwise velocity statistics for 
the coarse-grained case, which relies on results from ordinary 
CLPT presented in [18]. In the last subsection IV F we synop- 
tically compare the CLPT prediction for the streaming model 
ingredients to those of its coarse-grained generalizations scg¬ 
CLPT, cgCLPT and TCLPT. 


A. CLPT formalism 

In analogy to [18, 35] we define 


4i,/I 2 ) “ 



(38a) 


which allows to compute the ingredients of the streaming 
models arising from Z as given in (31) 


dPZ 

(idJi,)P 



ik-(a-r) f ddi dTz 

(27r)3 J (27r)2 


xF(Ai)F(A2)Kj,j,_,-^(k,q,Ai,d2). 


(38b) 


By integrating over A we obtain the bias parameters, which 
are expectation values of derivatives of the Lagrangian halo 
density held F[6x(q)] with respect to didq), according to [36] 


J ^F(A)(iAr exp , (38c) 

where cr^ - (6^(q)}. Furthermore also the integration over k 
can be performed analytically and only two dimensions of the 
q integration have to be done numerically, see [18, 35]. 

To numerically evaluate the real space correlation function 
(39a) and the pairwise velocity statistics (40a), (42a) within 
CLPT we resort to the C-H-i- code* that has been implemented 
by [18] for the dust model and which we extended to the 
coarse-grained dust case. The two quantities Kq and Ki do 
not depend on higher cumulants and hence do not dis¬ 

criminate between the two bias versions of cgCLPT (34/35). 
These quantities can be computed straightforwardly in full 
analogy to CLPT [18] by simply replacing dust correlators 
with their coarse-grained counterparts, which of course differ 
for TCLPT and cgCLPT. For K 2 , additionally the effect of ve¬ 
locity dispersion C® becomes relevant such that its form de¬ 
pends on whether scgCLPT (34) corresponding to s 0 

or cgCLPT (35) coiTesponding to is em¬ 

ployed. Both cases will be considered and compared to each 
other. 


• https://github.coi-n/wll745881210/CLPT GSRSD.git 
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B. Real space correlation function I + 


The real-space two-point correlation function 1 + ^x(t, t) is 
given by 




(27r)3 J (27r)2 

xF(Ai)F(A 2 )Ko(k,q,AuA 2 ), . 


(39a) 


where Kg has to be evaluated according to the cumulant ex¬ 
pansion theorem [53] 


Ko = \o(p 2 ) = exp 


Ati 

N=1 


(39b) 


0(Ph 


First we expand the exponent up to second order in the linear 
power spectrum Pp 


X At! (^^=o), 


v=o 


0(Pl) 

1 


- 1 - ^(^1 + ~ 2 ^ 12 ^ 2 ^L 


(Ai + A2)Uiki - -Aijkikj - -(d^ + d^)t/; 


20 ( 2 ) 


ki (39c) 


vll(2) 


ki - A 2 )A]fkikj - - Wijkkikjku, 


- iAiA 2 U; 

where the Lagrangian correlators are defined as 


jjmn(p) 

^mn(pqr) 


= = <^2)c , {6 i 62 )c 


= , A’;p‘’‘‘’ {6’^51 AY’A^J’)c , (39d) 


.mn{pq) _ 


:nxip)Aq)\ 




u 


and we adopt the shorthand notation f/, = t/*", A 
and Wijk = W™ introduced in [18, 35]. Whenever indices 
in brackets are omitted they have been summed over to the 
appropriate order of perturbation theory, for example f/, = 


-H t/® and An = h- A'^^^ 


+ ^(31)^ Note that 


‘J "ij ■ "ij ' "ij ■ “<7 

this notation for the correlators replaces that in Eq. (B4) which 
was defined according to [36]. In the following calculation we 
will keep this notation, such that whenever A, t/ or IT occur 
they refer to the usual kernels given in [18, 35]. In contrast, 
we will use A, t7 or IT for correlators arising from smoothed 
quantities 

= <52>. . h := { 6 i 62 )c , 1/'“® := (S'fdfAf), , 
{5'f6lAiAj\ , IT™ := . (39e) 

Plugging the exponential (39c) into the expression (39b) and 
keeping only the two terms exponentiated which are linear in 
the power spectrum and have non-zero limits as |q| —> cm gives 

-\(Fi+F,)d-\ 

1 

2 ' 


K{) = e 2 ‘J ‘ Je 


{l - AiA2^l - (^i + A2)Uiki + 


+ ^(^1 '^ifUikiU jkj + AiA2(Ai -I- A2)^LUiki (39f) 


{A\ + Aiwf^^hi - iAiAiU.'^^’ki 
(di H- A2)A]°kikj - ^Wijickikjkk] . 


Ai(2)- 


Note that this corresponds to Eq. (18) in [18] and resembles 
Eq. (72) in [35] except for the typo regarding the sign of the 
term. 

Eor the inference of the halo correlation function from the 
matter correlation function we used local Lagrangian bias (23) 
and fit the mass IgM determining the two bias parameters 
fei(lgMop() = {F'} and l72(lgMop() = <F") from (38c) that are 
listed in Tab. I for CLPT, cgCLPT and TCLPT. This procedure 
can be interpreted as fitting only bi while predicting b 2 ibi). 
We will give more details on the employed bias model and the 
fitting procedure in a forthcoming paper [46] . 


IgM 

13.00 

13.35 

13.59 

13.79 

13.99 

14.25 

14.67 

R(M)[Mpc/h] 

3.21 

4.20 

5.05 

5.89 

6.86 

8.38 

11.57 

CLPT (25) 

IgMop, 

12.93 

13.35 

13.63 

13.85 

14.07 

14.34 

14.79 

^lOgMpp.) 

-0.01 

0.26 

0.51 

0.77 

1.11 

1.68 

3.17 

^2(lgMpp,) 

-0.74 

-0.83 

-0.79 

-0.61 

-0.18 

1.06 

7.41 

cgCLPT (34/35) with or, = R(M) 

IgMpp, 

12.93 

13.35 

13.62 

13.84 

14.05 

14.30 

14.66 

i’lOgMpp,) 

0.00 

0.26 

0.50 

0.76 

1.08 

1.58 

2.65 

*2(lgMpp.) 

-0.74 

-0.83 

-0.79 

-0.62 

-0.23 

0.82 

4.70 

TCLPT (37) 

with 0 ^ = R(M) 

IgMpp, 

12.91 

13.33 

13.60 

13.83 

14.01 

14.26 

14.64 

^lOgMpp,) 

-0.02 

0.24 

0.48 

0.71 

1.01 

1.50 

2.58 

^2(lgMpp.) 

-0.74 

-0.83 

-0.80 

-0.66 

-0.32 

0.61 

4.35 


TABLE I: Best fit mass IgM^p, for the mass bins of average mass IgM 
for z = 0 with the corresponding bias parameters evaluated at IgM^p,. 


Pig. 4 shows the real space halo correlation function for all 
mass bins measured from the HR2 simulation in comparison 
to the linear and CLPT (25) predictions. In Pig. 5 we com¬ 
pare the smoothing effect on ^(r) caused by cgCLPT (34/35) 
corresponding to a coarse-graining in Eulerian space and our 
coarse-graining in Lagrangian space TCLPT (37) computed 
with a smoothed input power spectrum with a smoothing scale 
given by the Lagrangian size of the halo. That this smoothing 
scale should be relevant for halos will be discussed in [46]. 

As can be seen in Pig. 4, CLPT provides a quite accurate fit 
to the data points over a vast range of halo masses while only 
failing at small scales for the highest masses. As evident from 
Pig. 5 smoothing on the Lagrangian scale considerable flattens 
out the BAO peak around r llOMpc//; for both, cgCLPT 
and TCLPT, spoiling the agreement with the A-body data. It 
is therefore clear that keeping the smoothing scale and using 
a natural value for it, significantly affects the result even on 
scales that are naively much larger than the Alter size. 
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Plugging the exponential (39c) into the expression (40b) and 
expanding the second term gives 

X |!(/li + A2)Lfi + ikjAji - -(dj + Al)Uf'^ 

- AiA 2 Uj ' - ^-kjkk%i - (di + A 2 )kjA]f (40c) 

- idid2(di + A2)hUi - i(Ai + A2fkjUjUi 

- iAiA2^LkjAji - i{Ai + A2)kjki^UjA^j^ , 
where in addition to (39e) we defined 


FIG. 4: Normalized real space halo correlation function times 
for all 7 mass bins measured in HR2 (data points) and predicted from 
linear theory (thick black dashed) and CLPT (25) (thick solid) col¬ 
ored according to the mass bin. To achieve a clear representation for 
all masses we divided by the corresponding linear Eulerian bias 
hf(lgM^p, ,) and shifted all values by a constant 10(i - l)(Mpc//i)^ 
according to the i-th mass hin. 



350 


300 
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j. X ! CTq — /?(A/) 

—- cgCLPT: o-x = J?(M) 


r [Mpc//i] 


;= <5“5”A,.aA)- (40d) 

This result is analogous to Eq. (29) in [18] except for the use 
of coarse-grained instead of dust correlators. From this the 
pairwise velocity ui 2 (r) defined as r;i 2 ■ z - Vi 2 (r)r\^lr is com¬ 
puted according to (40). 

In Fig. 6 we show the CLPT (25) prediction for the pairwise 
velocity together with the HR2 data points normalized by the 
linear result, see [17], which is given by 




1 



dkPL(k)ii(kr). 


(41) 


In Fig. 6 we see that the CLPT prediction for the pairwise ve¬ 
locity is relatively inaccurate, especially compared to the ex¬ 
cellent agreement found in Fig. 4 for the real space correlation 
function. On the largest scales CLPT correctly describes the 
A^-body data and reproduces the linear theory result (41), how¬ 
ever around the BAO scale at 120Mpc//!, there is systematic 
offset. 


FIG. 5: Comparison between ^xf^ predictions of cgCLPT (34/35) 
(thick dashed) and TCLPT (37) (thin solid) smoothed on the La- 
grangian radius R(M) with the measurement from HR2 (data points) 
for the lowest 6 mass bins. For better visibility we normalized 
all functions by T'^'flgM^p,,) and shifted all values by a constant 
10(i - l)(Mpc//z)^ according to the i-th mass bin. 


C. Mean pairwise velocity vn 


The expression for the mean pairwise velocity v^ir, t) is 




(2;r)3 j (2;r)2 

X F(Ai)FiA 2 )Kij(k, q,AuA 2 ), (40a) 


where Ki is to be evaluated according to 
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FIG. 6: Mean pairwise velocity vi 2 from (40) compared to linear 
theory d’/j (41) for the 6 lowest mass bins for CLPT (25) (thin solid 
lines) and measurements from HR2 (thin dashed). 
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In Fig. 7 we show a comparison of cgCLPT (34/35) and 
TCLPT (37), both smoothed on the Lagrangian radius R(M), 
together with the HR2 measurements. 
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From Fig. 7 and Fig. 8 we deduce that cgCLPT systemat¬ 
ically increases the absolute value of the pairwise velocity 
compared to CLPT on small and on BAO scales. It shows only 
a small mass-dependence if the only the four smallest masses 
are considered. While cgCLPT improves the CLPT prediction 
on very large scales, TCLPT generally behaves better below 
90Mpc//!. TCLPT exhibits a stronger mass-dependence than 
cgCLPT. 


D. Mean pairwise velocity dispersion 


In order to evaluate the pairwise velocity dispersion cr 


^2 


12 


v\ 2 Vi 2 we have to determine 




12 


d/l 2 


(27r)2 

X F(Ai)F(A 2 )K 2 ,ijik, q, di, da), (42a) 


where K 2 is computed using Eq. (31). As mentioned before, 
K 2 depends on the tracer’s velocity dispersion which van¬ 
ishes identically for the models CLPT (25), scgCLPT (34) 
and TCLPT (37) but is relevant for cgCLPT (35). Therefore 
we split K 2 into one contribution K 2 from the dust model and 
another from velocity dispersion which is for cgCLPT 
(35) controlled by the smoothing scales oj and 07 , 


K 2 ,ij-K 2 ,j + K^^ 


(42b) 


The standard contribution to K 2 jj for the dustlike model is 
identical to Eq. (34-35) in [18] 


K 2 .ij = exp 


N=l ‘ JLiV=o ‘ 

“ :N+M 


N.M^O 


(43a) 
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FIG. 7: Comparison between predictions for the mean pairwise 
velocity P 12 compared to linear theory for cgCLPT (34/35) (thick 
dashed) and TCLPT (37) (thick solid) smoothed on the Lagrangian 
radius R(M) with the measurement from HR2 (thin dashed) for the 
lowest 2 mass bins. 
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FIG. 8: Detailed view on large scales of the mean pairwise velocity 
V 12 compared to linear theory for CLPT (25) (thin solid) and cgCLPT 
(34/35) (thick dashed) smoothed on the Lagrangian radius R(M) with 
the HR2 measurement (data points) for the lowest 4 masses. 


but evaluated with the kernels for the smoothed quantities 

K2ij — e^2'^'-'^'*.'e^2^2?+7;)D's 

X {(Ti + A2)^tJiDj - (Ai + A2)(A,ikkOj + Akjkktli) 

- AkikicAijki -H [1 - /liTg^L - (^i + A2)Uiikk]Aij 

+ i(Ai + A2)Ajj -t- iWkijkk^ , (43b) 

where in addition to (39e) and (40d) we defined 

A™ idl'd"AiAj}c , := (d"i'd"A^jAk),. (43c) 

In Fig. 9 the pairwise velocity dispersion cTj^ from CLPT 
(25) is shown together with the A-body measurements and 
the mass-independent result of linear theory, see [17], 

crl^ = 29{Y - ^2 dk Pdk) l^Mkr) - 

1 r“ 

J 

We can clearly see that CLPT significantly overestimates the 
amplitude of the pairwise velocity dispersion even on large 
scales and exhibits a strong mass dependence. The data points 
from HR2 lie close to each other for all masses which is 
captured by the linear theory result on large scales. The 
comparatively excellent performance of linear perturbations 
must be considered as accidental: the pairwise dispersion 
CTjj contains contributions of the one-point velocity variance 
{u'(x)d(x)), which we denoted by in the linear perturba¬ 
tion theory (44). This term is obtained as the limit r ^ 0 of 
{v‘(x + r)vdx)) and thus is sensitive to smallest scales. There¬ 
fore a large scale-independent offset error of cTjj is expected. 
In [17, 18] it has been accounted for that error by shifting the 
CLPT predictions to agree with the A-body measurements on 
large scales. 
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Fig. 10 compares scgCLPT (34) containing modified fluid 
dynamics to TCLPT (37) based on smoothing the input power 
spectrum. The smoothing on the Lagrangian scale R{M) 
significantly reduces the amplitude of cri 2 for all masses and 
narTows down the mass dependence compared to CLPT, see 
Fig. 9, bringing both models to better agreement with the data 
points. We therefore suggest that one should not adjust by 
hand the off-set of cri 2 , but should take this as an indication 
that a smoothing around the Lagrangian size of the halo 
should be applied. 

Contribution from the velocity dispersion of the tracer 
For the cgCLPT model (35) we have to consider the concep¬ 
tually new contribution to K 2 that arises from the velocity dis¬ 
persion of the tracers, encoded in according to 



FIG. 9: Pairwise velocity dispersion from (43) split according to 
(13c) into parts parallel cr^ and perpendicular a-\ to the pair separa¬ 
tion. Shown are the HR2 measurements (thin dashed) lying close to 
each other for all masses, the mass-independent linear theory predic¬ 
tion (44) (thick black) and the CLPT prediction (25) (thin solid). 



It is affected by both the spatial smoothing scale and the 
velocity smoothing scale 01 ,. Their contributions labeled 
and K'^‘ will be evaluated separately in the following. 

The corrections to the pairwise velocity dispersion cTj^ con¬ 
nected to the spatial coarse-graining are given by 



w=o ' 


J=0 


(1 -H 6)ViVj 
1+6 


- ViVi 



\(\+ 6)viVj 


\ 

+ 


(* 2 ( 92 )) 



[1+5 


L 


(xfq,)) (46) 


0(P]) 


and have to be calculated by explicitly evaluating the corre¬ 
sponding correlators as done in App. C. 

The correction term due to the coarse-graining cr,, with re¬ 
spect to velocity 


FIG. 10: Comparison between cTj^ predictions of scgCLPT (34) 
(thick dashed) and TCLPT (37) (thick solid) smoothed on the La¬ 
grangian radius R(M) with the measurement from HR2 (data points) 
for the lowest 6 mass bins. For better visibility we shifted all values 
by a constant 7(i - 1)( Mpc/h)^ according to the i-th mass bin. 
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can be obtained easily by combining (39c) and (39f) 

K‘^' = 2or,2 

X {l - 2AiA2h - 2(di + A2)Uiki - + Al)a-l 

- 2 ^ijkikj + -A\A^^\ + -did 2 (di + Al)^L&j( 

+ 3/li/l2(/li + A2)^LUiki + -A\A2^LAijkikj (47b) 

+ '^if'UikiUjkj + -(Ai + A2)(A^ + A\)&j(Uiki 

+ + '^ 2 )AijkikjUkkk - i(Af + 

- 2!di/l2t7"®^; - i(Ai + A2)A]fkikj - ^Wijkkikjkk]. 

In Fig. 11a comparison is shown between the predictions 
of scgCLPT (34), containing only modified dynamics corre¬ 
sponding to K 2 , and cgCLPT (35), also including velocity 
dispersion encoded in K^' and K^'. The inclusion of higher 
cumulants, which were assumed to be identical to those of 
dark matter, significantly increases the amplitude of the pair¬ 
wise velocity dispersion thereby spoiling the agreement with 
the data. This effect is due to the contribution of the spatial 
coarse-graining scale cTx while the velocity coarse-graining 01 , 
has hardly any effect for reasonable values of cr„. 



I I .. 

20 40 60 80 too 120 140 

r [Mpc/Zj] 


FIG. 11: Comparison between predictions from scgCLPT (34) 
{thin solid), based only on modified fluid dynamics encoded in K 2 
(43), and cgCLPT (35) (thick dashed), including higher tracer cu¬ 
mulants, smoothed on the Lagrangian radius R(M) together with the 
HR2 measurement (data points), upper panel The contribution 

(46) from the spatial smoothing o^. lower panel The contribution K^' 

(47) from the velocity smoothing oiJHq = 1 Mpc//i. 


E. Mean pairwise velocity skewness A12 


The first non-Gaussian correction A 12 in the ESM (12b) is, 
in analogy to 1112 and 0-^2 defined as 

[(. + &)A,«](...) = /<!=, / ^20- 

xF(Ai)F(A2)K3jjk(k,q,Ai,A2). 

(48a) 


As it was the case for K 2 , the explicit expression for Kj com¬ 
prises two parts 

K3,ijic = K3jjk + A 3 ";,.’I. (48b) 


The first K 3 jjk is the contribution from the dust model evalu¬ 
ated with the kernels for the smoothed quantities, while the 
second one only appears if higher cumulants of the 

tracer are present. To isolate the effect of non-Gaussianities in 
the pairwise velocity distribution we explicitly calculate only 
the dust part 


Ks.ijk = exp 
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FIG. 12: Pairwise velocity skewness A 12 from (48a) split according 
to (13d) into parts parallel A|| and perpendicular Aj^ to the pair sepa¬ 
ration vector. Shown are the measurements from HR2 {data points) 
which lie close to each other for all masses, the CLPT prediction (25) 
{thin solid) and the TCLPT prediction (37) {thick dashed) smoothed 
on the Lagrangian radius R = R{M). 

The CLPT result, when re-expressed in terms of the CLPT 
correlators defined in (39e), reads 

K 3 ,ijk = (faH)^ X 

X {iwuk + [(di ++ a<; v]} . (49) 

In Fig. 12 we show the effect of the pairwise skewness A 12 
which is the leading order non-Gaussian correction in the 
ESM (12b). Similarly as done for the pairwise velocity dis¬ 
persion cTjj, the skewness A 12 is also split into a part paral¬ 
lel and perpendicular to the line of sight according to (13d). 
The result shows that both CLPT and TCLPT based on the 
dust model fail to capture the non-Gaussian effects that are en¬ 
coded in the pairwise skewness basically for all scales below 
r - 100Mpc//;. However, it should be noticed that a smooth¬ 
ing on the Lagrangian scale does again bring the theory closer 
to the A^-body simulation, while reducing the spread in mass. 

F. Results 

The comparison of the streaming ingredients and the N- 
body data for large galaxy to cluster sized halos shows that 
the accuracy of the CLPT prediction is not sufficient to war- 
rent the use of the ESM, see Eig. 12. Hence, the GSM is suffi¬ 
cient when applied to CLPT based on the dust model and the 


ESM requires an improved theoretical modeling of the ESM 
ingredients with percent precision on 10 - 30Mpc//! scales. 

We investigated how the real space halo correlation func¬ 
tion and the velocity statistics are affected by a smoothing on 
the Lagrangian radius, the length scale corresponding to the 
physical size of a halo. In a forthcoming paper [46] we illus¬ 
trate that a smoothing on the Lagrangian size indeed optimizes 
the LPT prediction of the displacement field on the level of re¬ 
alizations compared to an A^-body simulation. Within our base 
model, GSM combined with CLPT, we observed a distinct be¬ 
havior depending on whether the coarse-graining is performed 
in the initial conditions (TCLPT (37)) or in Eulerian space 
resulting only in modified dynamics (scgCLPT (34)) or also 
including higher cumulants of the tracer (cgCLPT (35)), as¬ 
sumed to be identical to those of dark matter. Since all results 
have to be viewed in the light of fundamental limitations of 
perturbation theory to capture the correct small-scale behav¬ 
ior, we mainly focus on the impact of the smoothing on large 
scales for the model comparison. 

Fig. 5 illustrates the real space halo correlation function and 
that the dominant effect of coarse-graining on the Lagrangian 
radius is to damp the BAG feature around r a; llOMpc/A 
compared to CLPT. As expected, the impact of this modifica¬ 
tion becomes greater for massive halos since the Lagrangian 
size of the halo grows with mass. Given the already quite 
good agreement between the CLPT prediction and the data 
points evident from Fig. 4 and the strong damping, the La¬ 
grangian radius cannot be considered an optimal smoothing 
scale when aiming at optimizing the agreement with the N- 
body data for the real space correlation function. In the case 
considered in [46], where only the proto-halo centers in the 
realization are displaced, no damping around the BAG scale 
is observed. This means that the apparent disagreement most 
likely originates from the underlying local Lagrangian biasing 
scheme being not suitable to account for the fact that proto¬ 
halo centers are special rather than random points in the initial 
density field. We will address this issue in a future paper in 
the context of peak bias [54], where the correlation of density 
peaks is considered and spatial derivatives acting on the linear 
density correlation are included. We have preliminary results 
suggesting that the usage of peak bias indeed alleviates the is¬ 
sues present in the current model based on local Lagrangian 
bias. 

The pairwise velocity depicted in Fig. 7 shows that the two 
smoothing procedures either implemented in the initial condi¬ 
tions (TCLPT) or Eulerian space (cgCLPT) lead to different 
predictions. Smoothing only the initial conditions leaving the 
dynamics unaffected as done for TCLPT (37) increases the 
absolute value of the pairwise velocity compared to CLPT in 
agreement with the data. Smoothing in Eulerian space as done 
for cgCLPT (34/35) leads to a reduced mass-dependence and 
systematically increases the absolute value of the mean pair¬ 
wise velocity compared to CLPT, more prominently on small 
scales and around the BAG. cgCLPT improves the agreement 
with the data on large scales compared to CLPT for small and 
intermediate masses, see Fig. 8. 

For the pairwise velocity dispersion the main consequence 
of the coarse-graining on the Lagrangian scale is to downsize 
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the global amplitude and narrow down the mass dependence 
of the CLPT prediction as evident from Figs. 10 and 11. 
This effect is mainly caused by the mass dependence of the 
smoothing scale and there is only a minor difference between 
the modified dynamics contained in scgCLPT (34) and a 
smoothing of the initial conditions implemented by TCLPT 
(37). Furthermore, we showed in Fig. 11 that the cgCLPT 
model (35) including higher tracer cumulants signihcantly 
increases the amplitude of the velocity dispersion compared 
to the scgCLPT model (34) that follows the same modihed 
dynamics but does not include velocity dispersion. Most 
notably, the measurements from the A^-body simulation 
suggest that the bias assumption (30a) behind scgCLPT (34), 
neglecting higher cumulants of the tracers themselves, is 
more appropriate than cgCLPT (35) to describe dark matter 
halos. This is expected since proto-halos behave significantly 
more like a single-streaming fluid than the underlying dark 
matter distribution. 


V. CONCLUSION AND OUTLOOK 

In the first part of the paper we derived the Edgeworth 
streaming model (ESM (12)) for redshift space distortions 
which reduces to the Gaussian streaming model (GSM) when 
neglecting non-Gaussianities, starting from a general distri¬ 
bution function for dark matter tracers without making any 
assumption about the dynamics of those tracers and their re¬ 
lation to the dark matter field. We then studied the accuracy 
of both the GSM and ESM on the basis of the Horizon Run 
2 A^-body simulation halo catalog with masses ranging from 
large galaxy to cluster sized halos finding excellent agreement 
of the GSM for the quadrupole and hexadecapole on scales 
s > 30Mpc//! and an improvement of the leading order ESM 
over GSM on scales s < 30Mpc//r, see Eig. 3. 

In the second part of the paper we describe how to infer 
the streaming model ingredients for the GSM/ESM from per¬ 
turbation theory based on the dust model to determine the 
redshift space halo correlation function. Since halos carry 
an intrinsic length scale given by their Lagrangian radius, 
we considered a coarse-grained dust fluid to capture the ba¬ 
sic properties of the proto-halo fluid being comprised of ex¬ 
tended objects. In Sec. Ill we presented two different ways 
to incorporate a smoothing into the dust fluid description, im¬ 
plemented either in Eulerian or Lagrangian space, and pro¬ 
posed two possible biasing schemes to relate dark matter and 


halo phase space distribution functions. In Sec. IV we eval¬ 
uated the streaming model ingredients for the coarse-grained 
dust model within the Post-Zel’dovich approximation in its 
CLPT incarnation [35] which outperforms standard pertur¬ 
bation theory and even the Zel’dovich approximation in the 
nonlinear regime. We studied the impact of the smoothing 
on the streaming model ingredients and compared the coarse- 
graining in Eulerian space (scgCLPT (34) and cgCLPT (35)) 
to an approximate Lagrangian coarse-graining implemented 
by smoothing the initial power spectrum (TCLPT (37)). We 
found that halos can be described well as single-streaming 
tracers of the underlying dark matter distribution. Incorporat¬ 
ing a smoothing on the Lagrangian size of the halo turned out 
to have a strong effect on the BAG peak of the real space cor¬ 
relation function and hence as not suitable for improving the 
CLPT prediction when combined with local Lagrangian bias. 
By contrast, the smoothing improved the prediction for the 
pairwise velocity statistics, namely the mean pairwise veloc¬ 
ity for small masses and most notably the overall magnitude of 
the pairwise velocity dispersion irrespective of whether scg¬ 
CLPT or TCLPT are used. We expect that when combining 
the smoothing on the Lagrangian radius with the peak bias 
formalism [54] will considerably improve the agreement for 
the real space correlation function and the pairwise velocity 
statistics. 

In a forthcoming work [46] we investigate the appropriate 
choice of the smoothing scale to predict proto-halo displace¬ 
ment fields using the Zel’dovich approximation on the basis 
of A-body simulations. We then employ the fact that TCLPT 
can be implemented straightforwardly in CLPT by simply 
smoothing the input power spectrum and suggest a pragmatic 
optimizing scheme to accurately predict the ingredients of the 
GSM and hence the redshift space correlation function from 
the GSM whose individual accuracy on scales s > 30Mpc//! 
has been established. 
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Appendix A: Abbreviations 


abbr. full expression reference 

HR2 Horizon Run 2 halo catalog [38, 39] 

ZA Zel’dovich approximation [29] 

TZA truncated Zel’dovich approximation computed with smoothed input power spectrum [31] 

PZA Post-Zel’dovich approximation (higher order perturbation theory) [22] 

TPZA truncated Post-Zel’dovich approximation computed with smoothed input power spectrum [33, 34] 
LPT Lagrangian perturbation theory [22] 

CLPT Convolution Lagrangian perturbation theory (dust model), approximation to PZA (22), [35] 

cgCLPT coarse-grained Convolution Lagrangian perturbation theory (coarse-grained dust model) (34/35), [37] 
TCLPT truncated Convolution Lagrangian perturbation theory (dust model) (37) 

GSM Gaussian Streaming model (12a) 

ESM Edgeworth Streaming model (12b) 


Appendix B: Lagrangian correlators 


1. Lagrangian framework 


In Lagrangian perturbation theory the exact displacement field q) is expanded in a series with spatial parts and 

temporal coefficients given by the scale factor a{T) in an Einstein-de Sitter universe 

OO 

•P'(T,q-) = J]fl"(T)-l^<")(q-). (Bl) 

n=l 

The orders tfrC'O are expressed in Eourier space with the help of perturbative kernels in terms of the linear density field di 

= i J ^(27i^)3(«^P " ~ Pl...„)L^"^(Pl, . . ■,Pn)SLiPl) ■ ■ ■ SliPn) ■ (B2) 

Note that we employ here a different notation for compared to Eq. (A2) in [36] such that when translating the results an 
additional prefactor n! has to be taken into account. The vector valued kernels can be split into a longitudinal component 
and a transverse part according to 

LM ^ s(n) ^ j.(n) ^ ^ S^"\pi,. .. ,p„) = 0 , fc ■ T^"\pu ... ,p„) = 0 with k:=pi + ...+pn. (B3) 


In addition to those definitions, it is useful to define the following mixed polyspectra of the linear density field and the 
displacement field in the same way as done in [36] 


{SM ■ ■ ■ 6dki)Vt^\px) ■ ■ ■ ^ffiPm))^ 


{InfdYiiki + ... + ki +Pi + ...+ P;„)(-/)'"C}"/.7"'”^(fci ,...,ki;pi,...,pj, (B4) 


where an angle bracket with index c denote cumulants (connected correlators). Eor computations up to 1-loop level we only 
have to consider terms up to O(P^) which implies I + tii + n 2 +n^ <4 since due the properties of the cumulants only terms with 
I + m < 3 are relevant. Eurthermore only even I + ni + 112 + 11 ^ € 2N contribute because the initial density field is assumed to be 
a random Gaussian field. Eor / -H m = 2 we adopt the simplified notation 


C(k):^C(k,-k) , Ci(k) Ci(k;-k) , Cij(k) Cij(k,-k). (B5) 

The C as defined in (B4) should not be confused with cumulants of the phase space distribution function. We will also encounter 
mixed polyspectra with some '5' replaced by ’S' or by time derivatives. When combining the Gaussian streaming model with 
Convolution Lagrangian perturbation theory (CLPT) we will adopt another notation for the correlators, see Eq. (39d) in Sec. IV 
and [18, 35]. Eor the correlators involving time derivatives of ’S' we simply use that ’S'^"* = n/'TT’S'*"* and similarly for The 
specific index structure enforced by translation symmetry allows to describe Lagrangian correlators in real space entirely in terms 
of scalar functions of ^ = |q| = Iqi-qil, since any tensor can be decomposed in terms of and Eor example, any rank-1 tensor 
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can be written as Ti(q) - T{q)qi and similarly any rank;-2 tensor can be decomposed according to Tij(q) = Tsq{q)6ij + Tqq(q)qiqj. 
The ^-dependence of the functions can be expressed using spherical Bessel functions, namely 


x: 


d/i cos(xfi) - 2jo(x) 


x: 


dfi cos(xju) = 2 io{x) - 2 




We maintain the notation used in [36] and [35], in which the R and ^-functions, defined as 


Qnik) - —r I drPi(kr) I dx Piikyl - 2rx -i- r^) Qn{k, r, x ), 

Jo J-i 

f 


k^ 

4^2- 


Rn(k)^—Pdk) I drPL(kr)R„(k,r), 


(B6) 

(B7) 

(B8) 


have been computed for the standard fluid case. 


2. Lagrangian correlators for the coarse-grained dust model (cgCLPT) 

In the following we state the results for the Lagrangian correlators obtained from the coarse-grained dust model (35). Since 
for the CLPT computation, see Eqs. (B20-30) and (B41-46) in [35] only Qi, 2 , 5,8 are relevant, they are the only ones which will 
be listed here. 


Qi(k,r,x) 


Qi(k,r,x) 
Qiik, r, x) 
Qiik, r, x) 
Q 5 (k,r,x) 

Q 5 (k,r,x) 
Qiik, r, x) 
Qsik, r, x) 


-crJP(2j-2rx+l) 


36 (r2 - 2rx -H 1)^ 


49 1^1 + 3x^ + 4 (^-rx + r^x^ - 


lim 2i(^,r,x) = , 

^^-0 (r2-2 rx+lf 


- 14 [7(1 -H 3x^) 4 (-4rx -t 4r^x^ - lOrx^ -H 3r^x^)] 

-I- [ 49(1 + 3x^) -I- 36r^ + 4(-7rx - 1 Ir^x^ - 91rx^ + 51r^x^)j “''jj 

r2(x2-lf 


(B9a) 


(B9b) 


r(rx- l)e-oi"<:H 2 E- 2 «+i) 

6 (r2 — 2 rx + 1)^ 

rx(x^ - 1 ) (rx - 1) 


lim Q 2 {k, r, x) 

(77-^0 


[(7r - 2\rx^ + 34x^ - 20x) - (7r - 2\rx^ + 28x^ - 14x) (B9c) 

(B9d) 


re 


,crJk^(-Ap+2rx-2) 


6 (r2 — 2 rx -I- 1) 


(r2 — 2 rx -H 1) 

[(l3r -I- 14x + 34r^x^ — 20r^x — 41rx^) ^ — (ir + 14x + 28r^x^ — 14r^x — 35rx^) e”"' ^ 

(B9e) 


lim Qs{k, r, x) = 

(7[-^0 


(x^- l)(rx- 1) 
(r2 - 2rx + 1)^ 


^^-(Tjk'^(24~2rx+l) 

3 (r2 — 2rx + 1) 
lim Qs{k, r, x) = 


[(7rx2 - 7x) - (lOrx^ - 3r - 7x) 
r2 — 2rx + 1 


(B9f) 

(B9g) 

(B9h) 


Furthermore all quantities which contain the linear power spectrum, more precisely Eqs. (B20,25,41,44) in [35], have to be 
computed with the smoothed linear power spectrum, such that Piik) —» Prik) = exp Piik). 

In addition to the usual Rn=i, 2 , which are modified, we had to define another kernel Rq which accounts for the fact, that in our 
case the quantities cf'' and cannot be expressed in terms of ^ 1 . This is due to the different smoothing structure of cf^ and 
which both contain two quantities, compared to and which both contain three parts. As explained in [36] any 
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transverse part of is irrelevant such that we only obtain longitudinal parts for . for which, however, the transverse 

kernels have to be taken into account. 


This manifests itself in the following kernels which are instead of Eqs. (B26,43,46) in [35] then given by 


1 r°^ 5 

* 1 ■'* 2 l"“® 


2 Jiikq) 

3 kq 


1 r'^ 5 

= 5? 1 


-2jo(kq) + 6 


;i(M 


kq 


Ro(k, r) = 


480r3 


3 (r^ - l)' 51 (r^ + l) (eI - l)^] - Ei [-\o-,\\r + l)^]) 

r- 1 


•4 log 


r+ 1 


((7r2+2).-'^'=^^%2l(r2 + l)) 


168r(3/ - 5r4 + 9? - 3) - 8r(21r® - 50/ + 19? - 
1 


- 186^°^ ' 


(.(Txk? 
2 


(r + 1)(-17 + 51r + 44r^ + 68r" - 51r^ + 17r") 


+ 18e° 


{CTxkf 

8 

{(^xkf 

1 

{CTxk? 

2 

{CTxkY 


(17 + 44r + 88r^ + 68r" - \1?) 

192 

(17rH24r+ 11) + 


{o-xkf 


(r - 1)(17 + 51r - AA? + 68r' + 51r'‘ + 170 


(-17 + 44r - 88 ? + 68 ? + 11?) 


{(Txky 


■{\1? -2Ar+\\) + 


192 

((TxkY 


Ri(k, r) = 


-cr/k^(A+l) 


288? 


3{?-ir 


28 log 


r- 1 


r+ 1 


- i7eioi^*^('^+i) (Ei [-ko-??(r - l)^] - Ei [-ka??{r + 1)^]) 


-56r(-3rO llrO llr2-3) 
- 12 


- 12 


(o7k)2 

1 

(o-xk) 


(l7/ + 22? + 17) + [? + l) cosh [o??r) 


(BIO) 

(BID 

(B12) 


(B13a) 


(B13b) 


Ri(k, r) — lim R\{k, r) - lim R(){k, r) - 

(Tx-^i) fTv->0 


^ (r + 1) (17/ - 90r= + 17) - 0: (">■' + 22>-" + ") " O * '^1 } 

- 6 ? + 22 ? + 22 ? - 6r - 3 (D - 1)"^ log |^| 


48D 


(B13c) 


Riik, r) = 


(D - l)e-‘2iV(r2+l) 


288D 


3(r2-l) (rOl)p81og 
+ 56r(3D-2r0 3) 


r- 1 


r + 1 


(Ei [-\(Tx^k\r - 1)'] - Ei [-\ax^k\r + 1)"]) 

(B13d) 


- 12 

- 12 


34r 


(->) 


32r 

{(y-xkf 


cosh (cr??r^ 
34 


—^ ill? - 22? + 17 )-('■^ + 1 )-siiih (o??r) 1 

(o-xk)^^ > {a-xk?^ > {a-xk)(-\ V" 0 


R 2 {k, r) - lim R 2 (k, r) — 

(Tr^O 


(1 -?) 6 ? - A? + 6 r + 3[? - if [? + 1 ) log 1^1 


48D 


(B13e) 
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Note that, in the limit ^ 0 we correctly recover the result of [36]. 


Appendix C: Contribution to K 2 (42b) in cgCLPT 


In the following we explicitly state the explicit result for the oi-correction term to K 2 which is given by the second tracer 
cumulant 


Kh - exp 


Z 

N=\ 




y- 

^N\ 


N=0 


X 


j=o 


(1 + 5)viVj 
1+6 




(1 -H 6)ViVj 


ViVj 



( 3 ^ 2 ( 92 )) 



[1 + 5 




(Cl) 


Note that we used the notation y, for the dust velocity which yield the ordinary kernels for CLPT as given in [35], [18]. In 
contrast, we use v, (1 + 6 )vil(l + 6 ) for the mass-weighted velocity which corresponds to our cgCLPT kernels computed for 
v{x{q)) = a’S'iq). We have performed a similar calculation to the one presented as pedagogical example B4 in [35] and rely on 
results given in [36]. For convenience, the derivation was carried out by computing mixed correlators between coarse-grained 
density and displacements, 6 l and tP, and dust quantities, 6 l and S'. Since the extra terms encoded in Eqs.(Cl) are given in 
Eulerian space we first had to perform a mapping to Lagrangian coordinates. This has been done according to x{t) - q + ^(T,q) 
by using the Jacobian Fjj = dxijdqj - dq + Fij with determinant Jf - detE/j. More details concerning the mapping from 
Eulerian to Lagrangian space for the case of the coarse-grained dust model and explicit relations for the Lagrangian kernels up 
to third order can be found in [37]. The result up to 0(Pj) is 


k: 




where 


Kl 


r. 


2,ij 




-cr2((bl,2) 


- [exp (-a-y) [Ri(k) + 2 R 2 (k)] - [Ri(k) + 2 R 2 {k)]\ - 


1 , j likq) 

kq 


[exp [-cry] [Ri{k) + 2 R 2 {k)] - [R,{k) + 2 R 2 {k)]\ (yo(^?) - 3| §,§;•} 


tqfi 


r^(h0,2) f 

Jj - J 


kq 


Ak 

2 tF 


\Q. 2 (k) + £l 2 (k)] + - [Q 4 (k) + £l 4 (k)] + [Qi(k) + Qi(k)] Mkq) + [Q 2 (k) + D 2 (k)] ^ -j^j 


jiikq) \ 


(C3) 


(C4) 


dij (C5) 




^[Ei(k) -.Ro(^)] + ^[Qi(k) - Qiik)] 


6 ij + [Q 2 (k) + ^ 2 (k)]{jo(kq) - 3 


J ijkq) 

kq 


qtqj 


with the spherical Bessel functions jo and i\. The R and Q terms are the usual CLPT kernels from [36] whereas R and Q are our 
corresponding cgCLPT kernels given in (B9) and (B13). The additional kernels, Qi _4 and <3i_4, appearing in (C5) define two 
other classes of functions besides R and Q, according to 


n\ _ _ 

Qn{k) - J dr Piikr) J dx Pi(k Vl - 2 rx + r^) Q„(k, r, x) , (C6a) 

Q„(k,r,x) - —> 0 

£i„(k) -— 2 :rPi(k) f dr f dx Piik Vl - 2 rx + r^) Q„(k, r, x ), 

47r2 Jg 

~ (Tt—»0 

£i„(k,r,x) -—> 0, 


(C6b) 
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with 

Qi{k,r,x) - 
Qiik, r, x) - 
Qi{k,r,x) - 
Q/i{k, r, x) - 

&i(k, r, x) ■ 
& 2 (k, r, x) ■ 

^ 3 (k, r, x) ■ 
C!4(A:, r, x) ■ 




-cr//t^(2r^+l) 


Ik^ (r^ - 2rx + 1) 

[3rx^ -r-2x) 

Ik^ (r2 - 2rx + 1)^ 
r(l _ jc2^g-t5:V(2r"+l) 

Ik^ (r2 - 2rx + 1)^ 

[3rx^ -r-2x) e-°i^<=^(2'^+i) 


[(lOrx^ - 3r - lx) - lx (rx - 1) , 


(C7a) 

(C7b) 


Ik^ (r2 - 2rx + 1) 


[ 3 r (l - x^) e'^-^^^( 2'-^+0 + (lOrx^ - 3 r - 7 x) _ 7 ^ („ _ j) g 2 o;^/:Vxj ^ ^(. 7 ^^ 

[3r(l - + (lOrx^ - 3r - 7x)- 7x(rx - 1), (C7d) 


2r‘*(l _ ^2^ e-<^xV('^+2) 
7A:2 (r2 - 2rx + 1)^ 

2^3g-(7;^«:^(r^+2) 


[6 (x^ - 1) + (21rx - 34^^ + 6) - 7 (3rx - 4x^) , (C7e) 


[6(1- x2) (3rx2 -r-2x) 


7A:2 (r2 - 2rx + 1)^ 

+ (7r2 (5 - 9x^) X + 2r (5 Ix"^ - iV - 4 ) - 68x^ + 40x) 

+7 (r2 (9x2 _ 5 ^ ^ + 2r (-hx^ + x2 + l) + Sx^ - 4x) ^ 

2r4 (1 - x 2 ) e-‘^-'*'('^+ 2 ) 


(C7f) 


[6 (x^ - 1) e--^^^('-^+2) ^ ^21rx - 34x2 + 6) - 7(3rx - 4x2)e2‘"^^^^"] , (C7g) 


7A:2 (r2 - 2rx + 1)^ 

9„3 

[6 (1 - x2) (3rx2 - r - 2x) ^"-''='('•'+2) 


7^2 ^^2 _ 2 ^x + 1 )"^ 


+ (7r2 (5 - 9x2) + 2r (5 Ix^ - 19x2 - 4 ) - hSx^ + 40x) 

+7 (r2 (9x2 _ 5 ) ^ + 2r (-hx^ + x2 + l) + Sx^ - 4x) . 


(C7h) 



